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Abstract 

We demonstrate the viability of improved staggered light quarks in studies of heavy-light systems. Our 
method for constructing heavy-light operators exploits the close relation between naive and staggered 
fermions. The new approach is tested on quenched configurations using several staggered actions combined 
with nonrelativistic heavy quarks. Exploratory calculations of the Bg meson kinetic mass, the hyperfine and 
IP — IS' splittings in Bg, and the decay constant fs^ are presented and compared to previous quenched 
lattice studies. An important technical detail, Bayesian curve-fitting, is discussed at length. 
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I. INTRODUCTION 


Precise calculations of hadronic matrix elements are important ingredients in the quest to 
constrain the flavor-mixing parameters of the Standard Model, the Cabibbo-Kobayashi-Maskawa 
matrix elements Vff. For example, the main theoretical input in extracting the ratio |Vic;/hts| 
involves a combination of the decay constants, Jb and fs^, which parameterize leptonic B and 
Bg decays, and of the neutral B and Bg mixing parameters, Bb and Bb^- Uncertainties in these 
quantities, or more specifically in the combination ^ = (/s„ /fB)\/B bJBb, currently restrict our 
ability to carry out stringent consistency checks of the Standard Model (e.g. see [1]). If these 
theoretical errors could be reduced by a factor of 2 or more the impact would be immediate and 
far-reaching. Similarly, high precision theoretical calculations of form factors governing B D£i> 
and B Trii? decays are crucial to determinations of \ Vcb\ and \Vub\, respectively. 

Monte Carlo simulation of QCD on a lattice will ultimately provide the most accurate theoretical 
determinations of mixing parameters, decay constants, and form factors since lattice QCD is one 
of the few systematically improvable approaches to QCD. Understanding and removing systematic 
uncertainties in lattice calculations, however, is arduous and complicated, and much of the effort 
in lattice gauge theory over the past decade has focused on this task. One very promising outcome 
of all this activity is the emergence of improved staggered actions for light quarks combined with 
highly improved glue actions. The MILC collaboration, for instance, works with the “AsqTad” 
quark action which is free from the leading discretization errors, including those arising from the 
breaking of the fermion doubling symmetry, so that the action is accurate up to O (agO^) errors. 
They employ the one-loop Symanzik improved glue action with errors coming in only at O (a^a^). 
Staggered actions have an exact chiral symmetry at zero mass and are much cheaper to simulate 
than Wilson-type quark actions, so it has been possible for the MILC collaboration to carry out 
unquenched simulations with much smaller dynamical quark masses than has been attempted in 
the past. They are now starting to obtain impressive results for light hadron spectroscopy and 
light meson decay constants [2, 3]. 

In this article we demonstrate that improved staggered quarks can also be used very effectively 
to simulate the light quark in heavy-light systems such as in S-physics. The past decade has seen 
significant progress in our ability to simulate heavy quarks accurately (the commonly used NRQCD 
action, for instance, has errors coming in at O (agO^) and O (agA^cD/M), and work is underway to 
remove the latter). On the other hand, only Wilson-type actions have been used for the valence light 
quarks in heavy-light mesons, baryons, and electroweak currents, making it difficult to go much 
below mstrange/2 in the light quark mass due to the necessary computational expense. Consequently, 
the extrapolation of simulation results to the chiral limit leads to the dominant systematic error 
in studies of B and D mesons (aside from quenching uncertainties). Furthermore the leading 
discretization errors in heavy-light simulations come from the light quark sector since Wilson-type 
actions have worse hnite lattice spacing errors than improved glue, improved staggered, or NRQCD 
actions. This situation motivated us to initiate a new approach to heavy-light simulations, namely 
the use of improved staggered light quarks combined with nonrelativistic heavy quarks. Our 
approach can trivially be modified to use a Wilson-like action for the heavy quark instead of 
NRQCD. The goal is to simulate B physics at much smaller light quark masses than has been 
possible in the past and significantly reduce chiral extrapolation errors in decay constants, form 
factors and mixing parameters. Work toward this goal has already started on the MILC dynamical 
configurations [4]. It is important, however, to first establish that we understand how to combine 
staggered light and NRQCD or Wilson heavy fermions to form heavy-light operators, that we are 
able to carry out sophisticated fits to simulation data and extract physics reliably, and that these 
methods produce results in agreement with well-established results. It is for the last reason that 
this article focuses on the Bg system on quenched lattices, where methods existing in the literature 
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provide a solid basis for comparison. We present results for Bg meson kinetic masses, some level 
splittings, and the decay constant as evidence that our approach is working. For quick reference, 
we summarize the results of our finest isotropic lattice in Table I. Note that several systematic 
uncertainties remain, notably the error from determining the spacing on quenched lattices and 
discretization errors from using coarse lattices with the present level of improvement. Therefore, 
the results presented in this paper are useful for comparison to similar lattice calculations, but they 
are not appropriate for inclusion in phenomenological analyses. Having established this method 
as a promising approach, work is now underway on unquenched lattices to remove or reduce the 
systematic uncertainties and obtain state-of-the-art lattice QCD results. 

In the next section we introduce and describe the formalism for combining staggered light 
quarks with heavy quark fields to form bilinear operators that create heavy-light mesons or repre¬ 
sent heavy-light currents. A significant simplification comes about from recognizing the equivalence 
of staggered and naive fermions and writing down bilinears in terms of the latter. This will be 
explained thoroughly below. In Section III we give simulation details starting with a description 
of the glue, heavy quark, and light quark actions and then a discussion of our constrained fitting 
methods based on Bayesian statistics. Section IV gives results for heavy-light spectroscopy, in¬ 
cluding kinetic masses and a calculation of the Bg meson decay constant fss- Three appendices 
contain details regarding the theory, notation, and fitting techniques, respectively. 


II. FORMALISM 

In this section we describe how to combine naive/staggered light quarks with heavy quarks to 
form heavy-light meson and electroweak current operators. We adhere to the recently introduced 
practice of calling the doubler degrees of freedom “tastes” rather than “flavors” [5] (see also [6]). 
We will be guided by the following properties of naive/staggered actions. 

1. Up to overall normalization factors, there is no difference between using naive or staggered 
valence quarks in meson creation or current operators. Since naive fermions are easier to 
interpret and to handle theoretically, we will construct our heavy-light bilinears using naive 
fermion fields rather than staggered fields. 

2. Any correlator involving naive fermion propagators can be rewritten in terms of staggered 
propagators. Since staggered propagators are cheaper to calculate numerically, when it comes 
to actual simulations we will always work with expressions that have been converted to the 
staggered fermion language and involve only staggered (and heavy) quark propagators. 

3. The taste content of naive/staggered actions can be determined either in the coordinate or 
the momentum space basis. For heavy-light physics and for perturbation theory we find the 
momentum space interpretation to be more useful. 

We start by reviewing naive fermions and the identification of different tastes in momentum space. 
We will then explore the taste content of B mesons that appear when naive fermions are combined 
with heavy fermions. We assume that the heavy quark action has no doublers, as in NRQCD, or 
that doublers have been given masses of order the cutoff via a Wilson term, as in the Fermilab 
approach [7]. Heavy-light systems are much simpler than light-light systems since the heavy quark 
suppresses the taste-changing processes of the naive/staggered quark. 


3 



A. The free naive quark action 


Most of our discussion in this section will be for free unimproved naive fermions. Taste iden¬ 
tification and relevant symmetries survive the inclusion of gauge interactions and of the O (a^) 
improvement terms incorporated into the action that we actually use in our simulations (see next 
section for a description of the full action). The free unimproved naive fermion action is given by, 


with 



+ 


m 




T(x) 


V^T(x) = i [T(x + a^) - T(x - a^)]. 


( 1 ) 

( 2 ) 


We work with hermitian Euclidean 7 -matrices obeying { 7 ^, 71 /} = 26^^- It is well-known that the 
action (1) describes a theory with 16 tastes of Dirac fermions and that it has a set of discrete 
“doubling” symmetries, 


T(x) ^ e“'’"»MgT(x) 

■^(x) ^ e“-^®'^(x)Ml. (3) 

g is an element of G, the set of ordered lists of up to 4 indices, 


G = {g : g = ifii, g 2 , ■ ■ ■), < ^^2 < ■ ■ ■} ; (4) 

e.g. (2), (0, 3), and (0,1,2, 3) are elements of G, as is the empty set 0. The corners of the Brillouin 
zone are denoted by the 4-vector Tig such that 




TT 

a 


0 


e 5, 

otherwise. 


( 5 ) 


The Mg are transformation matrices 


M, = n (6) 

with 

(7) 

An illustrative way to reduce the taste degeneracy of the naive action is to diagonalize the 
action in spin space. Let $(x) and 4>(x) be a new set of 4-component spinor fields related to the 
original T(x) and T(x) fields via the Kawamoto-Smit [ 8 ] transformation. 

T(x) = D(x) $(x) T(x) = 4>(x) D(x)I ( 8 ) 


with 




fi=0 


( 9 ) 
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In terms of these new fields the naive fermion action takes on a spin-diagonal form, 


= d>(x) 


1 . 




m 


«>(x) } , 


( 10 ) 


where 


r/^(x) = (_i)(^o+...+a;^-i)/a 


( 11 ) 


Staggered fermions reduce the taste-degeneracy from 16-fold to 4-fold. The spin-diagonal form 
of Eq. (10) tells us it should be possible to do so, since each spin component of 4>(x) is independent 
of the other components. One way to proceed is to define 1-component fields xi^) through 

<h(x) = e(x)x(x). (12) 


The c-number spinor e(x) is usually chosen to be constant, and one ends up with the standard 
staggered fermion action for the fields xi^)- Reference [9] goes through a more rigorous and 
general method for reducing the number of independent tastes from 16 to 4 which does not rely on 
first going through the Kawamoto-Smit transformation. They exploit the symmetry (3) to place 
constraints among the 16 different tastes so that only 4 of them remain as independent degrees of 
freedom. (See also [10] which uses the Hamiltonian formalism .) 

Equations (8) and (10) allow us to derive the simple but important relation between the naive 
propagator Gq,{x]y) and the staggered propagator G^{x;y). One has, 

Eq. (8) ^ G^,{x;y) = H(x) G$(x;y) H(y)t (13) 

Eq. (10) ^ Gq,{x]y) = / 4 G^(x;y), (14) 

with /4 equal to a 4 x 4 identity matrix in Dirac space. This leads to 

G 4 ,(x;y) = D(x)D'f(y) x G^{x]y). (15) 


We use the identity (15) repeatedly in the present work to go from bilinears expressed in terms of 
naive fermion fields to correlators written in terms of staggered propagators. It can also be used to 
rederive familiar staggered correlators (e.g. for pions or rhos) starting from simple naive fermion 
bilinears. We emphasize that (15) is an exact relation even in the presence of gauge interactions; re¬ 
expressed as a relation between the inverse of the naive and staggered actions, respectively, for fixed 
gauge helds, it is valid configuration by configuration, and hence also for the fully interacting naive 
and staggered propagators. The relation (15) also holds for improved versions of naive/staggered 
actions. 

Before going on to discuss heavy-light bilinears, we end this subsection on basic naive fermion 
properties by reviewing the momentum space identification of naive fermion tastes. We continue 
to use the notation of [9]. The momentum space spinors are given by. 


X X 

with the inverse relation given by 


T(x) 


We use the notation, 



Jk,D JkeD (27r)"^ 


^(x)= / 

Jk,D 


f = [ 

A, Da AgD 0 (27r)^ 


(16) 


(17) 


(18) 
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where D denotes the full Brillouin zone, — ^ < /c/i < ^, and Dg just the central region, 
In terms of the momentum space spinors the free action (1) becomes 


/ i^ik) 

^i'yn-sm{kna) + m 

lk,D 

a 

L M J 


V’(fc) 


( 19 ) 


Using the 4-vectors Tig this can be written as. 


•So = J2f 'ip{k + Tig) 

g Jk,Df^ 


- M 


1 

— sm 
a 


{[k + 7Tg]ga) + m 


Ipik Tig) 


( 20 ) 


The next step is to define 16 new momentum space spinors q^{k) labeled by the elements g of the 
set G (4) 

q^{k) = Mg'il;{k + 7Tg) , q^{k) = Tp{k + 7rg)Ml; (21) 

the matrices Mg are those of ( 6 ). In terms of these new spinors, q^{k), and upon using the relation 

Mg'^gMl sin([A: + vr^J^a) = 7 ^ sin(/c^a), ( 22 ) 


the action Sq becomes 


5o = H / q^{k) 

g Jk,D(j 


sm{kga) + m 


I 


q3{k). 


(23) 


Eq. (23) clearly describes an action for 16 “tastes” of Dirac fermions. The sum over the 
elements of the set G can be interpreted as a sum over tastes. The doubling symmetry (3) which 
in momentum space becomes 


V'(/C) ^ Mgll}{k+Tlg) 

ip{k) 'ip{k + Trg)Ml, (24) 

takes one q^{k) taste into another up to possible sign factors, egi,g 2 = ± 1 , defined through 

^ 91^92 ~ ^ 91 , 92^9192 [ 9 ])- 

B. Heavy-light bilinears 

To discuss heavy-light bilinears we introduce heavy quark fields which can stand for either 
heavy Wilson or nonrelativistic fermions (for the latter case we will use the notation T//(x) ^ Q{x) 
in later sections with Q{x) a 4-component spinor with vanishing lower 2 components). The simplest 
interpolating operator one could write down for creating a B meson with a heavy quark field T//(x) 
and a naive antiquark field T(x) is 


Wb(x) = Th(x)75T(x) . 


(25) 


Let us analyze Wb{x) in 3-dimensional momentum space. To do so we introduce the 3D Fourier 
transformed fields 


v;(k,f) =a3^e-**^'^^(x,t) , ^(k,t) = a3^e**^'^^(x,t) 


(26) 
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and similarly for the heavy fields ^h- It is useful to introduce a subset Gg C G that involves only 
spatial indices /r ^ j = 1, 2,3. The full set G can be built up out of gg and gtgg with gg G Gg and 
gt corresponding to g = 0 (and Mg^ = iysyo)- In analogy with (18) we have 


lk,D, 


f d^k 
IkeD, (27r)3 


/k,D. 


f d?k 
Lgd,. (27r)3 


(27) 


where Dg denotes the full 3D Brillouin zone, — ^ < kj < and Dg 0 the central region, — < 

kj < Then, as is shown in detail in Appendix A, 


^>VB(x,t) = 


dko 




gsSGs J-7rf2a 27r 

{^^(k + 7rg^,t)75 \^Mlq3‘>{k,ko) + (-l)‘Ml^g^g®‘»'>(k,/ cq)] } . (28) 


For gg ^ the field V’_ff(k + creates a heavy quark with large spatial momentum so that 

any state containing it will have a large energy. Consequently, the contributions to the heavy-light 
bilinear from low-lying states come from the gg = 9 part of the sum in (28) 


r r/ 2 a 
J\<..Ds (^J—7T/2a 27r 



(k, t) 75 


g®(k,A:o) -h 



(29) 


In contrast, light-light bilinears receive contributions from all 8 sections of the spatial Brillouin 
zone (this can be seen by replacing by V’ in (28) and then using (21)). The gg ^ t) contributions 
to heavy-light bilinears are discussed in more detail in Appendix A, where we consider more general 
bilinears and show that they couple either to exactly degenerate states or to artificial high energy 
lattice states. 

Let us point out that in (29) there are contributions from both the pseudoscalar and the scalar 
state, which has a coefficient alternating in sign. The oscillating parity partner appears in light- 
light correlators as well. In Section III B we discuss how fits are able to separate these contributions 
from correlation functions. 


C. Heavy-light two-point correlators 

Once heavy-light bilinears with naive light quarks have been introduced, it is straightforward to 
obtain bilinear-bilinear two-point correlators and write them in terms of staggered propagators. 
Starting from this subsection we will revert to the usual practice of working with dimensionless 
spinor fields. Hence one should assume all 'L, h and y fields have been multiplied by a factor of 
and that all propagators are now dimensionless. Denoting the generic bilinear as 

Wr{x) (30) 


one has 


n'r„( 0 )) 


{r„ Gj, (0; i) G„{x-, 0)} 

X 

X c,c' 


(31) 
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where we have used eq.(15) to convert from Gq/ to G^. “Tr” stands for a trace over both 
color and spin indices, whereas “tr” stands for a trace over spin indices only. Using Gy.{0;x) = 

0)(—1)^^* one gets for the case F^c = ^sk = 75 

X 

= E E [tr (Sitx) G7 {i; 0)} G'/'=(x-, 0)] (32) 

X c,c' 

which couples to the B meson. For the B* meson, we set F^c = ^sk = 7 j which gives 

Cfi(p,t) = [tr{F!t(x)G^^(x;0)} G*/'^{x;0)] . (33) 

X c,c' 

In the above formulas we are now allowing the heavy-light mesons to have nontrivial momentum. 
As long as spatial momenta are restricted to apj < 7r/2 there should be no problems with the 
Lorentz and/or taste content of a meson suddenly changing at finite momenta. In later sections 
of this article we will present results showing good dispersion relations for B and B* mesons for 
momenta up to at least apj = tt/S to check this hypothesis. 

Although the discussion above implicitly assumes the use of local sources and sinks, generalizing 
to smeared sources and sinks is straightforward as long as one takes care that the smearing function 
preserves the doubling symmetry (3). This work employs local sources and sinks, with good 
results for the ground state mesons, but smearing is an important direction for future studies, 
especially those of excited states. Nonlocal sources have been used extensively in staggered fermion 
simulations of light hadrons. 


III. SIMULATION DETAILS 


A. Actions and parameters 


The gauge action used to generate the isotropic gauge configurations is the tadpole-improved 
tree-level O (a^)-improved action [11, 12] 

_ _ r, i ^ _ 1 R^ujx) _ 1 Ry^jx) \ , . 

I 3 12 12 I 

X,fJL>V K k ^ U fj. J 

represents the plaquette and the 2x1 rectangle in the (/U, v) plane; both are normalized 
so that {P^v) = {PfjLv) = 1 in the ^ oo limit. As part of our tests we also study anisotropic 
lattices where the temporal lattice spacing at is a few times smaller than the spatial lattice spacing 
Us] in this case improvement in the temporal direction is secondary to spatial improvement. The 
action used for the anisotropic lattices is the same in the spatial directions, but the rectangles 
with two units in the temporal direction are omitted and the space-time coefficients adjusted to 
be consistent with Symanzik improvement [13, 14]: 

^{aniso) _ ^ 1 / 


1 

lf>Pss'{x) 

1 (^) 

1 -^s^s (^) 'I 


13 ui 

12 

12 ul j 

&{ 

'^Pst{x) 

1 Rst{x)\ 


.3 ulul 

12 uiul J ■ 



( 35 ) 





For the values of the inverse coupling j3 and the bare anisotropy used in this work, the tadpole- 
improvement Landau-link factors Ug and ut, the spatial lattice spacing and the renormalized 
anisotropy ^ = ag/at were determined in Ref. [15]. The simulation parameters for the gauge 
configurations are summarized in Table II. 

The parameters for the isotropic lattices were intended to give approximately the same spatial 
lattice spacings as the anisotropic lattices. The isotropic 8^ x 20 lattice parameters were discussed in 
Ref. [16]. The isotropic 12^ x 32 configurations were generated for this work, and we determined the 
lattice spacing by calculating the static quark potential and using the phenomenological parameter 
ro = 0.5 fm [17] to set the scale. 

The light quark action we use is the O (a^) tadpole-improved staggered action [18, 19] which 
contains in place of the simple covariant difference operator in (1) an improved difference operator 
constructed as follows. First, the link matrices are replaced by “fat-link” matrices [20]: 

/ \ 

V,ix) = J] 1 + ^ U,{x) (36) 

\ ^ / symmetrized 

which contain 3, 5, and 7-hnk paths, all bent to fit within an elemental hypercube (Ref. [18] lists 
each term explicitly, and we write the second-derivative operator in Appendix B). This smear¬ 
ing effectively introduces a form factor in the quark-gluon vertex which suppresses the coupling 
of high momentum gluons to low momentum quarks. Second, the fat-link is further modified by 
adding what has come to be known as the Lepage term [19] in order to cancel the low momentum 
O {a?) error introduced by (36); 

F;(x) = iy(x) - ^ ^C/m(x) . (37) 

Finally, the remaining O (a^) (rotational) errors are subtracted by including a cube of the difference 
operator, the so-called Naik term [21]; therefore the O (a^) improved action is obtained by the 
replacement 

Vm ^ V;,-i(V^)3. (38) 

This action has been used in many recent simulations, quenched and unquenched, most prominently 
by the MILC Collaboration who call it the “AsqTad” action. In order to apply tadpole improvement 
consistently, powers of the covariant difference operators, (V^)" and {Vli are obtained by n 

/o'\ 

successive applications of or V)j , respectively, with no tadpole factors, replacing Ufj_/u^ 

in the final expression only after setting terms like U^{x)Uj^{x) equal to 1. In other words, one writes 
every operator in (36) in terms of paths of links, dividing each link variable by its corresponding 
tadpole factor u^. 

In this work we utilize anisotropic lattices, for which the improved staggered action is rewritten 
breaking the sum over spacetime directions into spatial and temporal parts 



The parameter cq is tuned to give the correct pion dispersion relation. We include a parameter 
yt,naik which we set equal to 1 or 0 whether we want to include the 3-link hopping in the temporal 
direction or not; we still call the yt,naik = 1 action “AsqTad”, and we refer to the yt, na.ik = 0 action 
as “AsqTad-tn.” Note that the isotropic AsqTad action is recovered by setting cq = ^ = yt, naik = 1 - 
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The NRQCD action is [22, 23] 


5NRQCD = E|<^tVt - 


Hq is the nonrelativistic kinetic energy operator, 


atHo = -- 


A( 2 ) 


2 ^{asMo) 

and 5H includes relativistic and finite-lattice-spacing corrections, 

^2 


at5H = -Cl ——-h C2 


8 ^{asMo) 8 (asMo) 

<^3 Tt — TTT2'^ • (V X E — E X V) 


• E - E • V 


8 {asMo) 

1 


-C4 ■ 


AW (A(2))2 

- ^ + Cb '77777~irr\ ~ *^6 


2^{asMo) ° 24^(asMo) ° 16n^2(asMo)^ 

All derivatives are tadpole improved and, 


f, A(4)=Evf) 

i=i i=i 


A( 2 ) = Evf, 


( 40 ) 


(41) 


(42) 


(43) 


Vfc = Vfc - iv® (44) 

The dimensionless Euclidean electric and magnetic fields are, 

2 ^ijkPij ■ (^^) 

Explicit expressions for m = 2,3,4 and are given in Appendix B. In most cases we set 

all 6 of the Cj = 1 and refer to this as the 1/M^ NRQCD action, even though the leading 1/M^ 
relativistic correction is also included. In order to make corresponding perturbative calculations 
simpler, some simulations were done setting ci = C 2 = C 3 = ce = 0 with C 4 = C 5 = 1 , and we call 
this the 1/M NRQCD action. In practice the results depend very little on which action is used, 
since the nonrelativistic approximation is very good for B mesons. 

The bare mass of the heavy quark, Mq, is chosen to be close to the bottom quark mass, based on 
simulations with Wilson-like light quarks [23, 24]. The bare mass of the staggered quark m is tuned 
to be close to the strange quark mass using the condition that the ratio of the “ss” pseudoscalar 
meson mass to the “ss” vector meson mass is approximately equal to ^J2m\ — = 0.673. On 

unquenched lattices the (/> mass is probably not accurately determined since it should be sensitive 
to the sea quark masses decreasing through the threshold for cj) ^ KK. Instead one should first 
determine the lattice spacing, then use 2m\ — m\ to determine the bare strange quark mass. On 
the other hand, for the quenched lattices in this work, the ratio serves as an appropriate fiducial 
for comparison between different lattices. 
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B. Fitting methods 


The light quark propagators are computed with anti-periodic boundary conditions in the tempo¬ 
ral direction; in contrast, the evolution of the heavy quark in time requires only an initial condition. 
Due to this difference, heavy-light meson correlators with temporal separations greater than Nt/2 
will be contaminated from the light quark propagating backward in time from the source across 
the time boundary, so we only compute heavy-light meson correlators up to Nt/2. The periodicity 
of the light quark can still be utilized to improve statistics by evolving the heavy quark backward 
in time from the source. We average the forward and backward propagating meson correlators 
configuration by configuration. 

The process of fitting the meson correlators to a series of exponentials is complicated because 
the temporal doubler causes the correlation function to couple not only to states with the quantum 
numbers expected from the continuum limit, but also to states with opposite parity times an 
oscillating factor (—Thus, we expect the meson correlators to have the form 

Kp — l Kp+Ko — l 

f{t-{Ak,Ek}) = £ + Y. (46) 

n=0 k=Kp 

which includes Kp states of expected parity and Kq states of opposite parity. In our study we always 
take Ko = Kp or Kp — 1 , and for the excited state energies we use the differences NEk = E^ — Ek-2 
as parameters in the fit. The K = Kp + Kq terms in the fitting function (46) can be rearranged as 

fit-{Ak,Ek}) = + (-l)*+Mie-^i* 

-h ^ ^^g-(ASfc-rA£;fc_ 2 -r...)i _ (- 47 ^ 

k=2 

Note that terms with even k are simple exponentials and those with odd k are oscillating expo¬ 
nentials. 

Recently a curve fitting method has been introduced to our community which allows one to 
estimate the systematic uncertainty from the series of states (47) in the correlator [25, 26]. One 
fits the correlation function C'(p, t) for all computed values of t, varying the number of terms K in 
the fit. For a given K, the best fit is obtained by minimizing an augmented x^- 

2K—1 \2 

xLs(CW,{A,},{(ft,i,))) = x"(C(t),{\}) + ^ ■ {«) 

i=0 

where we have generically denoted the parameters of (47) by 

A = (ylo, Eq, Ai,Ei,A2, AE 2 , ..., Ak, AEk) (49) 

the i-th component of which is Aj. 

In Appendix C we give a pedagogical summary of [25] as it applies to our calculation, but a 
few remarks here are in order. The second term on the right-hand side of (48) is the contribution 
of Gaussian priors for each fit parameter, and one sets the prior means and half-widths 5^ 
based on reasonable prior estimates for those quantities. The procedure is best illustrated by an 
example. Let us take a pseudoscalar heavy-light correlator, computed with the unimproved, or 
“1-link” staggered action, on the 8 ^ x 20 lattice as an example (see Fig. 1). The set of prior means 
/i and half-widths <5 used in fitting this correlator is given in Table III. The ground state energy 
and amplitude prior means were estimated from effective mass plots and the prior widths set at 
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50% and 25%, respectively. Priors for the excited states biased the amplitude fit parameters to be 
all of the same order and the energy differences to be equal and about 300 MeV, roughly the size 
of the 2S — IS* and IP — 15* splittings in B spectrum. Recall that the NRQCD action does not 
include the rest mass, so the energy illsim is equal to the physical meson mass minus an energy shift 
A. Tables IV-VI show the results of fits to the propagator in Fig. 1 as the number of exponentials 
changes from 2 to 8. The uncertainties are estimated from the inverse of the matrix VVXaug of 
second derivatives ([VVy^uglii = 


= 



(50) 


which assumes the shape of near its minimum (A* = Af“, Vj) is quadratic in A* 


Xaug Xauglinin 


5 


min' 

i y 




^ Aaug 

dXidXi 


{Xj-Xfn 


(51) 


In Fig. 2 we plot the non-oscillating and oscillating ground state energies, as well as the hrst 
excited state energy, vs. the number K of exponentials in the fit. The rest of the fit parameters are 
given in Tables IV and V. One can clearly see the stability of the ground state fit parameters Aq, Eq 
and Ai,Ei as K is increased. The beginning of a plateau at iF = 3 implies at least one excited 
state is needed in the fit in order for the excited state effects to be removed from the ground states. 
Table VI similarly indicates that iF > 3 is necessary in order to have an “acceptable” x^^g/DoF; 
as we discuss in Appendix C, y^ug/DoF should only be used as a gross check of the fit. E.g. 
Xaug/DoF > 2 implies the fit function is a highly improbable model of the data, but one should 
not necessarily prefer a fit with Xaug/I^oF = 0.8 over one with Xaug/I^oF = 1.3. 

Note that the uncertainties estimated from the fit for the ground state parameters are much 
smaller than the widths of the corresponding priors, while the errors from the fit are comparable 
to the prior widths for most of the excited state parameters. The hrst excited non-oscillating state, 
/c = 2, is an exception, appearing to be well constrained by the data until another non-oscillating 
state, /c = 4, is included in the ht. This means that the iF = 3 and iF = 4 ht result for E 2 does 
a good job of absorbing the effects of the excited states, but that there is not enough constraint 
from the data (or the priors) to separate the hrst excited state from the second. Thus, we conclude 
that iF = 3 is sufficient to obtain reliable estimates of the ground state energies and amplitudes 
and that the data is not sufficiently precise to extract excited state energies and amplitudes. 

We are able to utilize this constrained curve htting method to ht all of our data except in one 
case: the heavy-light correlators computed with the AsqTad action on the 8^ x 20 lattice. We were 
not able to hnd hts with Xaug/^^o^ V 8; one example is shown in Fig. 3 where the ht is visibly 
much worse than for the l-link action shown in Fig. 1. This turns out to be a consequence of using 
an action with next-to-next-to-nearest neighbor couplings in the t-direction on a lattice with coarse 
temporal lattice spacing. 

The free Naik fermion dispersion relation (see Fig. 4) has complex solutions which implies 
there may be excited states with negative norms contaminating the correlators at short time 
separations. If the temporal extent of the lattice were sufficiently long and sufficiently precise 
correlators were computed, these negative norm states which have energies proportional to 1/a 
would have a negligible effect; one could include only points with t greater than some minimum 
value in the fit, or one could include a negative norm state in the fit. However, for the 8^ x 20 
lattice where 1/a = 0.8 GeV we are unable to drop enough points and get a good fit while keeping 
enough to fit to states of both parities. Also, when we tried to include a negative norm exponential 
in the fit, large cancellations with the positive norm excited states resulted in unstable fit results. 


12 



We checked this hypothesis on the 1/a = 0.8 GeV lattice by simulating with 4 different staggered 
quark actions: the 1-link and improved actions as well as an action where the Naik term was 
included but no fattening of the links was done (the Naik action) and an action with fat-links but no 
Naik term (the “fat-link” action). We were able to obtain reasonable fits to heavy-light correlators 
with the 1-link and fat-link actions, but not with the Naik action nor the fully improved action; we 
tabulate typical values for y^^g/DoF in Table VII. Furthermore we performed simulations on an 
anisotropic 8^ x 48 lattice with a very fine temporal lattice spacing l/a* = 3.7 GeV using the 1-link 
action, the AsqTad-tn action {yt,na.ik = 0), and the full AsqTad action {yt,na.ik = !)• In all 3 cases 
we found acceptable fits with similar values of Xaug/I^ol^j again tabulated in Table VII. Fig. 5 
shows the pseudoscalar propagator on this lattice for the AsqTad action. We have no problem 
fitting to heavy-light correlators on a finer isotropic 12^ x 32 lattice where 1/a = 1.0 GeV with the 
O (a^) improved action. Therefore, the origin of the problematic fits on the 1/a = 0.8 GeV lattice 
is due to a particular lattice artifact arising from the temporal Naik term; but with a larger lattice 
scale 1/a > 1.0 GeV these artifacts become insignificant. 

Let us return to the subject of estimating the uncertainties of the fit parameters. The second 
derivative of Xaug (^0) gives a reliable estimate of the uncertainty assuming that the priors are 
reasonable and that the data are approximately Gaussian. Resampling methods, such as the 
jackknife or the bootstrap, can be used to check whether the distributions are Gaussian, and they 
provide a simple check on statistical correlations between different fit parameters. Both procedures 
take many subsets of the data as estimates of the original set; performing a fit on each subset yields 
a distribution for each fit parameter from which an error estimate can be made. We employ the 
bootstrap method of resampling which requires some modification in order to properly handle the 
contributions of the priors: as we show in Appendix C one must randomly select new prior means 
for each bootstrap fit [25]. Table VIII shows the results of applying this bootstrap analysis to 
the heavy-light pseudoscalar correlator computed with the 1-link staggered action on the 8^ x 20 
lattice. These results can be compared to those in Tables IV and V. We find both methods 
produce comparable error estimates. For ease of error propagation, we use bootstrap method to 
quote uncertainties in the results presented below. 


IV. RESULTS 

This section contains several results produced using the methods proposed and described above. 
The purpose of this study was to check the validity of this proposal, so the results presented below 
should not be construed as state-of-the-art calculations to be used for phenomenology. The results 
here show that NRQGD-staggered calculations produce results comparable to NRQCD-Wilson 
calculations - central values agree and statistical and fitting uncertainties are comparable - but at 
a fraction of the computational cost. A more complete calculation of the B spectrum and decay 
constant on finer, unquenched configurations is underway which will exploit the advantages of 
improved staggered fermions to produce, we believe, the most accurate theoretical computation of 
those quantities to date. 


A. Light hadron masses and dispersion relations 

As mentioned before we chose a value for the bare staggered mass m so that the ratio of the 
light pseudoscalar mass to the light vector mass would be somewhat near the phenomenological 
value 2m\- — = 0.673. We use the pseudo Goldstone boson correlator Tr|G^(x;0)p to 

compute the pseudoscalar meson mass and the correlator (—1)*'=/“ Tr|G^(x; 0)p to compute the 
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vector meson mass. These masses and their ratio are listed in Table IX for the different lattices 
and actions. Note that even on the a~^= 0.8 GeV lattice the light hadron correlators from the 
AsqTad action do not suffer the contamination from the negative norm states which affected the 
heavy-light correlators, as discussed in Sec. IIIB. 

One measure of discretization effects is the dispersion relation. Specifically, we can compute 
the “speed-of-light” factor 


“ ln|2 


(52) 


which should equal 1 in the absence of lattice artifacts. The Naik term (38) is responsible for 
subtracting the O [a?) uncertainties in c^(p) and its success can be seen in the following results. 
Table X lists the values of (? computed with several values of momentum (averaged over all equiv¬ 
alent orientations in momentum space). On the coarser lattice {(3 = 1.719) one can see that using 
fat links does not improve (? much compared to the 1-link action, but adding only the Naik term 
to the 1-link action results in a significant improvement. This is borne out on the finer lattice 
{j3 = 2.131), where the AsqTad action has an improved (?. Fig. 6 shows comparison of c^(p) 
between these results and those for improved Wilson actions [16]. The AsqTad action has a better 
pion dispersion relation than the clover action, but not quite as good as the D234 action. 

On the anisotropic lattices we use this quantity to tune the bare parameter cq in (39); it is 
adjusted so that the pion speed-of-light parameter c^(p) ~ 1. Table XI lists the values of cq we 
used and the resulting computed values of (? for several momenta. 


B. Finite momentum Bs and its mass 

The energies, £'sim(p)) extracted from correlation functions include an unknown but momentum 
independent shift due to the neglect of the heavy quark rest mass, i.e. 

Flsim(p) = ^(p) - A (53) 

where E(p) is the physical energy. In perturbation theory, the shift A is the difference between 
the renormalized pole mass and the constant part of the heavy quark self-energy; 

Apert = ZmMo — Eq . (54) 

Given Fisijn(O) from a simulation, the physical mass of a hadron can be computed through 

Mpert = Esimi^) + ZmMq — Eq (55) 

where we attach the label “pert” to denote that the perturbative shift Apert was used. For the 
finer isotropic lattice and the 1/M NRQGD action with aMp = 5.0 we find 

Z]\/[Mq — Eq = Mq — 0.890 (Xg + Mq O (^6) 

The results obtained on the finer isotropic lattice, using the AsqTad staggered action, give 
Mpert(7?s) = 5.55 ± 0.45 GeV and Mpert(7?*) = 5.58 ± 0.45 GeV. The numerical size of the 
O (a^) uncertainty is estimated by taking ag ~ 0.3, a typical value for quenched lattices with 
these spacings, and assuming the coefficient of the O (a^) term is 1 (times Mq as indicated in 
(56)). 

The physical mass can also be calculated nonperturbatively, using the dispersion relation 

F;2(p) = + |p|2. (57) 
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In order to cancel the unknown shift in (53), we consider (-E(p) — -E(O))^ 
which we square and solve for the mass 

_ IpP ~ [-E'sim(p) “ -£^8101(0)]^ 

" 2[i?sim(p)-i^sim(0)] • 


(£;sim(p) - -Esim(O))^, 


(58) 


When the mass is computed using (58), we call it the kinetic mass, to distinguish it from the 
perturbative result Mpert- Setting |p| = 27r/12a = 0.52 GeV, the kinetic masses on the hner 
isotropic lattice (with the AsqTad light quark action) are = 5.56 ± 0.33 GeV and 

Mkin(.B*) = 5.68 ± 0.33 GeV. 

Figs. 7 and 8 show the kinetic masses for the Bg and B^ for several momenta. We hnd excel¬ 
lent agreement between the perturbative and nonperturbative calculations of the mass. Further¬ 
more, the consistency of the kinetic masses over several momenta demonstrate that the combined 
NRQGD-improved staggered formulation gives the correct dispersion relation for Bg and B^ up to 
IpI = 7 r/3a = 1.1 GeV. One should not put too much weight on any agreement or disagreement 
between the calculation and experiment, given that the calculation is quenched, the lattice spacing 
not precisely determined, and the quark masses not precisely tuned. 


C. Mass splittings in the Bg system 

Since the shift A between simulation energy and the physical energy (Eq. 53) is entirely due 
to the NRQCD action, it is universal for all bound states with the heavy quark. Therefore, we 
can compute mass splittings much more precisely than suggested by the uncertainties in M^in and 
Mpert- The splittings we compute on various lattices which correspond to the Bg system are given 
in Table XIV; below are a few remarks concerning the different calculations. 

The hyperfine splitting M^j — Mb^ is the most straightforward to compute since it is the 
difference between the Esim for the non-oscillating ground states of the vector and pseudoscalar 
correlators. The results are comparable to previous quenched lattice studies; Figure 9 shows our 
quenched results on the 2 isotropic lattices compared to results published in Refs. [27, 28, 29, 30]. 
This splitting was also computed using NRQCD in Ref. [31], but they have different systematic 
errors caused by the quenched approximation, specifically they set the bare bottom quark mass by 
tuning the T mass, instead of a heavy-light mass. Our error bars are larger than those for most 
other results for two reasons. The first is simply that this work is based on 200 configurations 
compared to 300 [28], 278 and 212 [29], and 2000 [30] (Ref. [27] used 102 configurations). The 
second is that the Bayesian curve fitting method includes as part of the quoted uncertainty an 
estimation of the error due to excited state contamination, in contrast to the single exponential 
hts used in previous work. 

Quenched results have an inherent ambiguity depending on which physical quantities are used to 
set the lattice spacing and bare quark masses. Preliminary results on unquenched lattices indicate 
that the inclusion of sea quarks yields a unique scale and bottom quark masses [32] and give a 
Bg — Bg splitting [4] consistent with the experimental measurement Mb* — Mb, = 47.0 ± 2.6 MeV 
[33]. 

The L = 1, or “P-wave”, states B*q and Bgi have the same quantum numbers as the oscillating 
states in the pseudoscalar and vector correlators, respectively. The fact that Egim for these states 
can be computed using the same correlator data as the L = 0 states should be another advantage 
over formulations with Wilson-like quarks. In practice, however, it appears that the coupling of 
these states to the local-local correlator is rather small and consequently the fitting uncertainties 
for these splittings are large. Smeared sources and sinks for both heavy and light quark propagators 
should be explored as methods for amplifying the coupling to the P-wave states. In Table XIV we 
list some combinations of splittings. 
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D. Decay constant 


The heavy-light decay constants are defined through the matrix element of the electroweak axial 
vector current 


(0|To|S^) = (0| 975706 | 5 s) = IbsMbs - 


(59) 


The fields in the current above are those defined in the Standard Model, so a matching must be 
performed between them and the fields of our lattice action. The continuum heavy quark field 6 is 
related to the nonrelativistic field cj) through the Foldy-Wouthuysen-Tani transformation 



where 


Q = 



(61) 


Expanding the QCD axial-vector current in terms of NRQCD operators up to O (Agco/dT) and 
at O (as) in perturbation theory yields a combination of the three operators 

= ^7570 g (62) 

-^Q (63) 

Jo^ = 7570 <3- (64) 

The operator equation is then written as 

^0 = (1 + (^sPo)Jo ^ + (1 + a^pi)Jg ^ -|- asP2JQ ^ ■ (65) 


The symbol = is meant to imply that matrix elements of the operators on the left and right hand 
sides are equal, up to whatever order in the effective theory we are working. Since we are neglecting 
terms of order asAgco/Af the terms proportional to pi and p 2 are dropped from our analysis. The 
relation we use to do the matching is 


Ao = (1 + «.po) 


( 66 ) 


where the 1 /(aM) power law mixing of Jg^^ with Jg^^ is absorbed at one-loop level into a subtracted 
Aqcd/M current [34] 


T(l,sub) 


4 "^ - asCio 


(67) 


and po - Cio = Po- 

Since the heavy spinor obeys 7570(3 = the matrix element (0|Jg°^|i?s) is related simply to 

(2i) 

the ground state amplitude of the pseudoscalar heavy-light correlator {p = 0, t). Let us denote 
this amplitude by Coo, then 


Coo 




( 68 ) 
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To get the Kqqd/M current matrix element we compute correlators where we put at the sink. 
Let us denote the ground state amplitude of this correlator by Cio, then 


Cio = 






2Mb^ 


(69) 


As mentioned before, we concentrate on the quenched 12^ x 32 lattice which is closest to the 
target unquenched configurations, albeit coarser. Fits to these correlators, shown in Fig. 10, yield 
the bootstrapped ratio 




- 0.064 ±0.005 (stat). 


(70) 


The 1/M NRQCD action is used for this calculation, for which we compute (with aMo = 5.0) 
Pq = 0.208 ± 0.003 and Cio = —0.0997. Performing the subtraction (67) we find 


±l,sub)| „ , 

^ ° = -0.034 ±0.004 (stat). (71) 

{0\Jr\Bs) 


This ratio can be compared to other lattice formulations; it is the “physical” Aqcb/M correction 
to with 1/a power law effect subtracted at the one-loop level. The 3.4(4)% correction we 
find on the a~^ = 1.0 GeV lattice is in excellent agreement with the 3-5% corrections found using 
the NRQCD and clover actions on lattices with inverse spacings from 1.1 - 2.6 GeV [34]. Note 
that even on the finest lattice in Ref. [34], where power law contributions are the largest, the 
one-loop subtraction takes (Oj jRs)/(0j = —14% to (Oj jRs)/(0l jRs) = —4%, 

in agreement with calculations on coarser lattices. Given the present agreement between our result 
and that of Ref. [34], we can expect a similarly successful subtraction in our ongoing calculation 
with the unquenched MILC ensemble. 

Applying (66) and (59) gives the quenched result 


fBs = 225 ± 9(stat) ± 20(p.t.) ± 27(disc.) MeV. (72) 

The 20 MeV perturbative uncertainty is the estimate of the O (a^) error in (66), obtained by 
taking a* 0.3 and a coefficient equal to 1. The other perturbative uncertainties, due to one-loop 
corrections to the coefficients in the action and in the operator matching, are O (asAgcD/Af) which 
is estimated to be 2.4%, assuming Aqcd = 400 MeV (and using Mq = 5.0 GeV). The 27 MeV 
discretization error is our estimate of the ©(a^aAgco) error in the current (62); again we 
assume the leading correction term comes with a coefficient of order 1. This error may be reduced 
to O (Q;s(aAQCD)^) by improving a la Symanzik, which requires calculation of pi and p 2 in 
(65) [35, 36]. Finally, note that we have included the O (ag/aM) power-law correction; we would 
have estimated this to be a 6% effect, but it was calculated to be 3% (compare 70 and 71). Given 
those uncertainties, we find agreement with the recent quenched world average Jb^ = 200 ± 20 
MeV [37]. 


V. CONCLUSIONS 

We believe the methods outlined within this paper provide the quickest route to accurate cal¬ 
culations of B meson masses and decay constants on realistic unquenched lattices. Improved 
staggered fermions have several advantages over Wilson-like fermions and are far less expensive 
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to simulate than domain wall or overlap fermions. The equivalence between naive and staggered 
fermions greatly simplifies the construction of operators which couple to states of interest. The fact 
that the NRQCD action does not have a doubling symmetry leads to taste-changing suppression 
in heavy-light mesons, avoiding the ambiguities of the light staggered hadrons. 

We have presented results on several types of lattices, the most important being the finer of 
the two isotropic lattices since it is most similar to the unquenched MILC lattices. The results 
from these simulations have no unpleasant surprises: they agree with results produced by previous 
quenched simulations. Therefore, we can trust this formulation when it is used in parts of parameter 
space inaccessible to other formulations. 
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APPENDIX A: FORMALISM DETAILS 

In this Appendix we present a more detailed analysis of the heavy-light operators used in the 
numerical calculation. 

Since naive fermions have 16 taste degrees-of-freedom, there is the possibility of forming 16 
different B mesons labeled by the light taste index g, i.e. Bg. The general choice for a Bg meson 
interpolating heavy-light operator takes on the form 

WBg{x) = ^H{xh5Mge"^x^^{x) (Al) 

The 16 different operators lead to degenerate states, since they are related by the symmetry 
transformation (3). It is sufficient to work with just one of the 16 choices to extract all the relevant 
physics. In our simulations we usually use the simplest choice 5 = 0, i.e. Eq. (25). Any other choice 
would have served equally well. For instance, consider the case g = gj with gj equal to one of the 
spatial directions and carry out a sum over spatial sites. Eq. (Al) then becomes 

^Wbj{x) = ^ (A2) 

XX X 

One sees that the zero spatial momentum Bj meson operator is identical to an operator one would 
superficially (and incorrectly) associate with a B* meson of polarization “j” with momentum tt/o 
in the direction. The correct interpretation of (A2) is that it represents a zero momentum 
pseudoscalar heavy-light meson. This will become more evident when we look at the operator 
y^Bg{x) in momentum space. We have verified that the RHS of (A2) gives identical correlators, 
configuration-by-configuration, to (25) (the latter summed over space). (In fact, the symmetries of 
(3) provide excellent tests of one’s simulation codes.) Therefore, it is sufficient to work with just 
one type of B meson operator, e.g. with just (25). 

In order to delve further into the Lorentz quantum number and taste content of the interpolating 
operators yVBg{x) we will look at this operator in momentum space. In terms of the “tilde” fields 
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(26) one has (we take the case where g does not include a temporal component; the latter case can 
be discussed in a completely analogous way) 


= / Vi/f(k,t)75Mg^(k + 7rg,t) 


= E 

Qs^Gs 




Vi/f(k + TTg^, t) ysMg ^(k + TTg + TTg^, t) 


(A3) 


We extract the taste content of this bilinear by writing 

dko ,, 

^ —Trja 


^{k + 7Vg + 7Vg^,t) = j ^e*^o*V'(k + 7rg + 7rg^,A:o) 

J —IV/a 


,ikot 


dko 
I —Tvlla 27r 

I—7vl2a 27r 


V'(k + TTg + 7rg^, fco) + (-l)V(k + TTg + 7Vg^,ko+ 7rgt) 
Mlgq<^^3^k,ko) + {-iyMlg^gq<^^>^‘S{k,ko) 


(A4) 


so that 


n/2a 

“^0 ^ikgt 


= E / / 9 


{^^(k + TTg^, t) ysMg /cq) + ko)] } (A5) 


Since there is no doubling symmetry for the heavy quark action, the field V’H(k + 7 Tg^,t), for 
TTg^ 7 ^ 770 , represents a heavy quark with large spatial momentum. Consequently, even though 
the operator in (A5) couples to zero momentum meson states, the states corresponding to g ^ 9 
are very energetic. This is precisely the important difference between studying heavy-light and 
light-light mesons with light staggered quarks. 

We will estimate the effect of the 7 ^ 0 sectors below, however the lowest energy state, and 
consequently the dominant contributions to a Wsgi^, t) correlator, will come from the region gs = 9 
in the sum • 


E (x, t) 

X 


dko i 


ikot 


k,Ds 07r/2a 27r 

{^H(k,t )75 q^(k,ko) + {-lYMlq^*<^{k,ko) } 


(A6) 


The non-oscillating contribution is the Bg meson of taste g. Its parity partner is a O"*" meson, 
usually called the J = 0 P—state. It is remarkable that both S— and P—states can be obtained 
from a single correlator. Note also that the combination ■0j:^(k, t) 75 (/^(k,/ cq), with its obviously 
pseudoscalar Lorentz structure holds for all tastes g, i.e. for trivial and nontrivial Mg in (Al). 

We discuss next those terms omitted upon going from (A5) to (A 6 ). Take for instance the 
contribution from gg ^ gl = gi with “r’ equal to one of the spatial directions. The non-oscillatory 
term becomes 


+ '^9ht)l5Mgiq^'-^{k,ko) = ^j|^(k-F tt^;, t) iq; g^'^(k, fco). (A7) 


One sees that the Lorentz structure is that of a 1“ particle. However, since the heavy quark has 
very high momentum and no doublers, this intermediate state is highly virtual. Such states would 
appear in fits to correlation functions as extra structure at energies of order AE ~ These 
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lattice artifacts can also affect low energy states through loops; their effects can be estimated 
perturbatively and are part of the O [asO?) errors inherent in the action. Such errors can be 
removed, if need be, by perturbatively improving the action further, but there is little evidence 
that they are important at practical values of the lattice spacing. 


APPENDIX B: DISCRETE DERIVATIVES AND FIELD STRENGTHS 


Here we write explicitly the higher order tadpole-improved derivatives and improved field 
strength tensor used in the fermion actions. 


The 


= —[U^{x)'^{x + a^) + Ul{x - a^)'^{x - a^)] 

- 2T(x) 

Vjf)T(x) = ^^[Uf,{x)U^{x + a^)'if{x + 2a^) 

- U^ix - a^)Ul{x - 2a/,)T(x - 20^,)] 

- — [Uf,{x)'i/{x + a^) - C/^(x - o^) T(x - a^)] 

V^'^^T(x) = ■\[Uf,{x)U^{x + a^)^{x+ 2 a^) 

+ Ul{x - a^)Ul{x - 2a^)^{x - 2a^)] 

-A — [U^{x)^{x + a^) + C/^(x - a^) T(x - a^)] 

Ufj, 

+ 6T(x) 

covariant derivatives acting on link matrices are defined as follows: 


u 


^ V^C/j.(x) = Uf^{x)U,y{x + af_,)Ul^{x + a,. 




Tv® UA^) = 

UjIJ 


- ul{x - a^) Uu{x - a^) U^{x -a^ + a^) 
^U^{x)Uu{x + a^)Ul{x + ay) 


uf^Uy L 


f7^(x - a^) Uy{x - a^) U^{x -a^ + Oy) 


-f7,y(x) 

Uy 


(Bl) 


(B2) 


(B3) 


(B4) 


(B5) 


The field strength operator F^y{x) is constructed from the so-called clover operator Vt.^y{x) 
Fii,y[x) = — (Q,^y{x) — Qj^i^(x)'j , 


t^uKX) - 

^fiu{x) = .22 Ua{x)Ui3{x + aa)U-a{x + aa + af3)U-i3{x + af3) (B6) 

where the sum is over {(a, P)}^y = {(/U, v), {v, —fj,), (—^, — 12 ), (— 0 )} for ^ / z/ and f7_^(x-|-a/i) = 
Ulix). The O (a^) improved field strength tensor is 

1 


Ff,y{x) = ^Ff,y{x) - ^ ■^{U^{x)Ff,y{x + a^)Ul{x) 

L 
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+ Ul{x - a^)F^^{x - a^)U^{x - a^)) - {fi ^ u) 

+ ^ “I 2 “ 

0 ur 


ut 


(B7) 


APPENDIX C: FITTING DETAILS 

In this Appendix we give a pedagogical discussion of the constrained curve fitting proposed in 

[ 25 ]. 

Recall that the standard fitting procedure is to minimize the function, or equivalently, to 
maximize the likelihood of the data, C{t), given a set of fit parameters. The likelihood probability 
is given, up to a normalization constant, by 


'P{C{t)\f{t;X),I) oc exp(-y 


where I represents any unstated assumptions. Explicitly, 

X' = - f{f,\)) K-] [{C{t')) - . 


(Cl) 


(C2) 


The correlation matrix, K, is constructed to take into account correlations between C{t) and C{t'): 

{{c{t) - {C{t))) (C(t0 - (C(tO))) . (C3) 

with N equal to the number of measurements. 

Usually one cannot include enough terms in the fit to account for excited state contributions 
before the algorithm for minimizing breaks down. The minimization algorithm diverges as 
it searches in directions of parameter space which are unconstrained by the data. In the past 
the solution has been to limit the number of fit terms, then discard data by including C{t) for 
t > ^min > 0; the optimal value of fmin is selected by a combination of looking for y^/DoF = 1, 
maximizing the confidence level (Q factor), and observing plateaux in effective masses. A major 
weakness of this procedure is that it provides no estimate of the error due to omitting the excited 
states from the fit. 

The constrained curve fitting method of [25], by using Bayesian ideas, allows one to incorporate 
the uncertainties due to poorly constrained states by relaxing the assumption that there are only 
a few states which saturate the correlation function. Bayesian fits maximize the probability that 
the fit function describes the given data, written as V{f{t;\)\C{t),I); this probability is related 
to the likelihood (Cl) by Bayes’ theorem 


Vif{t-,X)\C{t),I) = V{f{t;X)\I) 


ViC{t)\f{f,X),I) 

V{C{t)\I) 


(C4) 


and is called the posterior probability. The denominator in (C4) is treated as a normalization 
and plays no role in finding an optimal set of fit parameters. On the other hand, the prefactor, 
V{f{t; X)\I), which multiplies the likelihood is the prior probability; its inclusion is what permits 
fits to many parameters. 

The prior probability contains whatever assumptions about the values of the fit parameters one 
can safely make without looking at the data. In our case of fitting meson correlators, before any 
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fitting is done one has an idea of a range of possible values for the amplitudes and energies 
Ej^. Given such a range, the least informative prior distribution is a Gaussian with mean fi and 
half-width 5, in which case the prior probability is given by 

(C5) 

We sometimes refer to the set of as the “priors.” The quantity which is minimized in the 

hts is 

2K-1 Q _ x2 

xLg(C(*). {>.}.{(». •Si))) = ,y(C(t),{Aj)) + (C6) 

i=0 

oc -21nlP(/(t;A)|C(t),I) . 

Expression (G6) highlights how the prior distribution stabilizes the minimization algorithm. As 
one increases the number of fit parameters the terms from the prior in (G6) give curvature to Xaug 
which prevents the minimization algorithm from spending much time exploring the flat directions 
of in order to find a minimum. The trick now is to distinguish which parameters are constrained 
by the data and which ones are fixed by the priors. 

A remark on counting the net degrees-of-freedom (DoF) of the fit is in order. As usual each 
of the data points represents one DoF, but then each parameter of the fit which is constrained 
by the data uses up one of those degrees. However, in the Bayesian curve fitting method, there 
are several ht parameters which are unconstrained by the data and do not count against the net 
DoF. Usually there are a few parameters obviously constrained by the data and a few obviously 
determined solely by the prior, but there may be some parameters for which such a distinction is 
not clear. Therefore, we simply take the DoF to be the number of data points; instead of striving 
for a fit which produces x^/DoF < 1, we look for Xaug/I^oF 1 together with the property that 
the ratio stays constant as more ht terms are added. 

Given a sample of N measurements, one bootstrap sample is obtained by selecting N measure¬ 
ments, allowing repetitions, from the original N measurements. In principle one would perform 
a ht on every possible bootstrap sample generating a Gaussian distribution of bootstrapped ht 
parameters {A^}, the half-width of which gives the bootstrap uncertainties However, there 
are a total of {2N — 1)\/{N\{N — 1)!) ~ ways to make a bootstrap sample,^ so it is impossible 
to generate the entire bootstrap ensemble - it also unnecessary. The bootstrap distribution can 
be reliably estimated by randomly generating Nb bootstrap samples for large enough Nb- We use 
Nb ^ N = 200, and have check that changing Nb by a factor of 2 makes no signihcant difference 

In the unconstrained htting method x^ would be minimized for each bootstrap sample, resulting 
in a set of ht parameters which reproduce the likelihood probability distribution (Gl). For the 
constrained hts where we minimize Xaug) however, it is not enough to resample the likelihood - 
one must resample the whole posterior distribution, i.e. the product of the likelihood and the prior 
distribution. Therefore, for each bootstrap sample we randomly choose a new set of prior means 


^ Counting the number of possible bootstraps is equivalent to counting the number of ways n indistinguishable balls 
can be distributed into k distinguishable buckets: each bucket represents an original measurement and the number 
of balls in a bucket indicates the number of times the measurement appears in a given bootstrap sample (and 
n = k = N). The answer is called the integer composition of n into k parts and is equal to the binomial coefficient 
/ n + k — 1 \ 

V y 
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{;U^,}, b G [IjiVe] using the same distribution used in (C6) 


nA) 


1 



A - 





(C7) 


The bootstrap fits then yield an ensemble of results for each fit parameter with a nearly 
Gaussian shape 


r{\,) ~ exp(^-j <C8) 

where is the fit result for Aj on the 6-th bootstrap sample, (Ai)s = and a^. is the 

bootstrap uncertainty for A*. In practice one finds that the distribution of fit results is Gaussian 
shaped in the center but has stretched tails which artificially inflate the quantity ^ (Aj)^ — (A ?)b 
making it a poor estimate of cr^, [38, 39]. Instead we estimate the width of the bootstrap distribu¬ 
tion by discarding the highest 16% and lowest 16% of A^ and quoting the range of values for the 
remaining 68% as 2(T^,. Having obtained bootstrap fits to several correlators, say {A^} and 
we estimate the uncertainty in functions of the fit parameters g{Xi, Uj), for example mass ratios, by 
computing the function for each bootstrap sample and truncating the resulting distribution just 
as discussed for individual fit parameters. 
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TABLE I: Summary of quenched results from the isotropic 12^ x 32 lattice (1/a = 1.0 GeV). Results are 
checks of our new formulation, not state-of-the-art computations to be used for phenomenology. 


quantity 

result 


5.56 ±0.33 GeV 

MpeABs) 

5.51 ±0.45 GeV 

Mki„(R:) 

5.68 ±0.54 GeV 

MpeABt) 

5.53 ±0.45 GeV 

B* — Bg splitting 

25.0 ±4.8 MeV 

fs. 

225 ± 9(stat) ± 20(p.t.) ± 27(disc.) MeV 


TABLE 11: Simulation parameters for the quenched gauge configurations. There are 200 configurations for 
each parameter set. 


volume 

P 

A 

1/a, (GeV) 

as/at 

Us 

Ut 

UsMo 

8^ X 20 

1.719 

- 

0.8 

1 

0.797 

0.797 

6.5 

8^ X 48 

1.8 

6.0 

0.7 

5.3 

0.721 

0.992 

7.0 

123 X 32 

2.131 

- 

1.0 

1 

0.836 

0.836 

5.0 

123 ^ 4g 

2.4 

3.0 

1.2 

2.71 

0.7868 

0.9771 

4.0 


TABLE 111: Gaussian prior means /r and widths S for fits to pseudoscalar heavy-light propagator on the 
8^ X 20 lattice, m = 0.18. We use ek to denote Ek for the ground states (k = 0,1) and AEk for the excited 
states (fc > 2). 


k 

MAfc ± Ak 

Met ± 5et. 

0 

0.94 ±0.47 

0.900 ±0.225 

1 

0.94 ±0.47 

1.40 ±0.35 

2 

0.60 ±0.60 

0.40 ±0.30 

3 

0.60 ±0.60 

0.40 ±0.30 

4 

0.60 ±0.60 

0.40 ±0.30 

5 

0.60 ±0.60 

0.40 ±0.30 

6 

0.60 ±0.60 

0.40 ±0.30 

6 

0.60 ±0.60 

0.40 ±0.30 
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TABLE IV: Dependence of fit results on number of terms (K) included in fit - energies of the 8^ x 20 
heavy-light pseudoscalar correlator. Uncertainties quoted here were estimated from Wy^ug described 
in the text. 


K 

aEo 

Non-oscillating terms 
aAE2 

aAE^ 

gAEq 

2 

1.044(0.003) 




3 

0.919(0.013) 

0.492(0.048) 



4 

0.921(0.013) 

0.505(0.061) 



5 

0.915(0.018) 

0.365(0.161) 

0.282(0.250) 


6 

0.917(0.018) 

0.372(0.165) 

0.304(0.259) 


7 

0.914(0.022) 

0.322(0.199) 

0.269(0.257) 

0.348(0.297) 

8 

0.915(0.021) 

0.328(0.203) 

0.288(0.261) 

0.361(0.297) 



Oscillating terms 



K 

aEi 

aAUg 

aAE^ 

aAE’j 

2 

1.290(0.015) 




3 

1.503(0.028) 




4 

1.470(0.099) 

0.405(0.296) 



5 

1.461(0.103) 

0.422(0.293) 



6 

1.461(0.100) 

0.412(0.299) 

0.412(0.299) 


7 

1.455(0.102) 

0.419(0.299) 

0.419(0.298) 


8 

1.464(0.098) 

0.419(0.300) 

0.418(0.300) 

0.418(0.300) 


TABLE V: Dependence of fit results on number of terms (K) included in fit - amplitudes of the 8^ x 20 
heavy-light pseudoscalar correlator. Uncertainties quoted here were estimated from VVy^ug described 
in the text. 


K 

^0 

Non-oscillating terms 
A 2 

A 4 


2 

1.955(0.010) 




3 

1.047(0.104) 

1.088(0.093) 



4 

1.063(0.111) 

1.082(0.091) 



5 

0.999(0.174) 

0.599(0.411) 

0.557(0.411) 


6 

1.014(0.176) 

0.601(0.406) 

0.557(0.409) 


7 

0.980(0.224) 

0.518(0.437) 

0.505(0.456) 

0.185(0.445) 

8 

0.993(0.225) 

0.527(0.432) 

0.488(0.457) 

0.205(0.452) 



Oscillating terms 



K 

Ai 

^3 

^5 

Ay 

2 

0.478(0.013) 




3 

0.674(0.024) 




4 

0.580(0.263) 

0.105(0.289) 



5 

0.553(0.265) 

0.140(0.293) 



6 

0.583(0.268) 

0.008(0.448) 

0.119(0.324) 


7 

0.572(0.267) 

-0.005(0.448) 

0.159(0.336) 


8 

0.611(0.275) 

-0.043(0.460) 

0.005(0.480) 

0.179(0.412) 
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TABLE VI: Augmented chi-squared per degree-of-freedom for the fits in the preceding 2 tables. 


K 

xLr/DoF 

2 

47.7 

3 

0.60 

4 

0.69 

5 

0.66 

6 

0.74 

7 

0.83 

8 

0.92 


TABLE VII: Summary of fits to pseudoscalar heavy-light correlators. (“AsqTad” implies yt,naik = 1 unless 
otherwise indicated.) 


P 

action 

1/a, (GeV) 

l/ot (GeV) 

K 

xL./DoF 

£;sim(p = 0) (MeV) 

1.719 

1-link 

0.8 

0.8 

3 

0.59 

735(10) 

1.719 

AsqTad 

0.8 

0.8 

3 

8.93 

- 

1.719 

Naik 

0.8 

0.8 

3 

17.6 

- 

1.719 

Fat-link 

0.8 

0.8 

3 

0.51 

691(20) 

1.8 

AsqTad 

0.7 

3.7 

4 

1.59 

790(36) 

1.8 

AsqTad (tmin = 3) 

0.7 

3.7 

4 

0.87 

791(39) 

1.8 

AsqTad-tn 

0.7 

3.7 

4 

1.03 

901(19) 

2.131 

1-link 

1.0 

1.0 

4 

0.48 

873(9) 

2.131 

AsqTad 

1.0 

1.0 

4 

0.96 

765(9) 


TABLE VIII: Bootstrap ht results for the 8^ x 20 heavy-light pseudoscalar correlator for fits to K terms. 


A 

K = 2, 

K = 4: 

K = 5 

^0 

1.043(0.116) 

1.061(0.116) 

1.003(0.183) 

aEo 

0.919(0.014) 

0.920(0.014) 

0.918(0.019) 

Ai 

0.680(0.025) 

0.559(0.270) 

0.538(0.262) 

aEi 

1.508(0.030) 

1.457(0.113) 

1.450(0.112) 

^2 

1.098(0.102) 

1.098(0.101) 

0.736(0.329) 

aAE2 

0.499(0.057) 

0.514(0.070) 

0.405(0.127) 

^3 


0.141(0.313) 

0.178(0.297) 

aAE^ 


0.442(0.288) 

0.447(0.305) 




0.496(0.421) 

aAE^ 



0.380(0.247) 
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TABLE IX: Light pseudoscalar and vector meson mass, computed with the same light quark propagators 
used for heavy-light mesons. (“AsqTad” implies j/qnaik = 1 unless otherwise indicated.) For comparison, we 
nominally associate the physical strange sector with mps/mv = 0.673. 


p 

l/as (GeV) 

action 

atm 

mps (MeV) 

my (MeV) 

mps/mv 

1.8 

0.7 

1-link 

0.04 

843(6) 

1251(11) 

0.674(0.005) 

1.8 

0.7 

AsqTad 

0.04 

626(18) 

989(31) 

0.632(0.022) 

1.8 

0.7 

AsqTad-tn 

0.04 

628(19) 

994(37) 

0.630(0.024) 

1.719 

0.8 

1-link 

0.18 

761(1) 

1171(30) 

0.649(0.017) 

2.131 

1.0 

1-link 

0.12 

825(2) 

1218(27) 

0.678(0.015) 

2.131 

1.0 

AsqTad 

0.10 

685(3) 

1035(23) 

0.662(0.016) 

2.4 

1.2 

1-link 

0.03 

808(6) 

1108(21) 

0.728(0.012) 

2.4 

1.2 

AsqTad-tn 

0.03 

619(8) 

913(17) 

0.679(0.013) 


TABLE X: Speed-of-light parameter squared for several values of a|p| on the isotropic lattices. Since no 
tuning is done, one can estimate the size of lattice artifacts in finite momentum states from how different 
is from 1. 

P 

L 

action 

c^(27r/L) 

c2(2^/27^/L) 

(? j L) 

1.719 

8 

1-link 

0.656(6) 

0.631(8) 

- 

1.719 

8 

fat-link 

0.729(17) 

0.712(20) 

0.684(24) 

1.719 

8 

Naik 

0.883(12) 

0.916(22) 

- 

1.719 

8 

AsqTad 

0.892(14) 

0.922(22) 

- 

2.131 

12 

1-link 

0.794(8) 

0.767(26) 

0.778(16) 

2.131 

12 

AsqTad 

0.946(14) 

0.952(26) 

0.836(80) 


TABLE XI: Speed-of-light parameter squared for several values of a p 
parameter cq in the anisotropic action (39) is tuned so that « 1. 

on the anisotropic lattices. The bare 

P 

L 

action 

Co 

L) 

<P{2^/2t:/L) 

c^{2\/Z'k jL) 

1.8 

8 

1-link 

1.1 

1.004(23) 

0.993(31) 

0.991(46) 

1.8 

8 

AsqTad 

1.4 

0.940(46) 

0.952(56) 

0.975(48) 

1.8 

8 

AsqTad-tn 

1.4 

0.948(54) 

0.957(51) 

0.980(45) 

2.4 

12 

1-link 

1.0 

0.965(35) 

0.957(51) 

0.937(60) 

2.4 

12 

AsqTad-tn 

1.0 

0.931(44) 

0.957(51) 

0.853(81) 


TABLE XII: Bootstrap fit results for the 12^ x 32 heavy-light pseudoscalar {Bg) correlator for several 
momenta. (AsqTad light quark action.) 


A 

a|p| = 0 

op = 27r/12 

a|p| = 2^/27^/12 

a|p| = 2v^7r/12 

op = 47r/12 

^0 

0.129(0.009) 

0.130(0.010) 

0.129(0.011) 

0.126(0.012) 

0.139(0.017) 

citEo 

0.774(0.008) 

0.799(0.009) 

0.822(0.010) 

0.845(0.011) 

0.872(0.014) 

Ai 

0.071(0.045) 

0.071(0.045) 

0.071(0.046) 

0.070(0.050) 

0.070(0.053) 

at El 

1.298(0.090) 

1.312(0.095) 

1.326(0.095) 

1.343(0.100) 

1.370(0.089) 

A 2 

1.209(0.014) 

1.214(0.014) 

1.222(0.015) 

1.232(0.016) 

1.218(0.019) 

Ot Aii/2 

0.636(0.010) 

0.613(0.010) 

0.589(0.010) 

0.565(0.010) 

0.540(0.010) 

^3 

0.738(0.046) 

0.745(0.047) 

0.752(0.048) 

0.760(0.051) 

0.760(0.052) 

atAEs 

0.149(0.127) 

0.141(0.133) 

0.135(0.120) 

0.124(0.123) 

0.095(0.114) 
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TABLE XIII: Bootstrap fit results for the 12^ x 32 heavy-light vector (B*) correlator for several momenta. 
(AsqTad light quark action.) 


A 

a|p| = 0 

a|p| = 27r/12 

a|p| = 2v^7r/12 

a|p| = 2v^7r/12 

a p = 47r/12 

Ao 

0.112(0.011) 

0.112(0.013) 

0.109(0.016) 

0.104(0.018) 

0.122(0.022) 

atEo 

0.799(0.010) 

0.823(0.011) 

0.846(0.013) 

0.868(0.016) 

0.898(0.019) 

Ai 

0.070(0.046) 

0.070(0.048) 

0.069(0.049) 

0.070(0.052) 

0.069(0.052) 

at El 

1.339(0.103) 

1.348(0.083) 

1.364(0.084) 

1.379(0.093) 

1.394(0.079) 

^2 

1.232(0.014) 

1.238(0.015) 

1.247(0.015) 

1.258(0.018) 

1.240(0.022) 

atAE2 

0.599(0.009) 

0.575(0.009) 

0.552(0.010) 

0.530(0.011) 

0.503(0.010) 

^3 

0.747(0.047) 

0.754(0.049) 

0.761(0.051) 

0.767(0.052) 

0.767(0.052) 

atAE^ 

0.116(0.121) 

0.113(0.110) 

0.105(0.107) 

0.096(0.101) 

0.078(0.092) 


TABLE XIV: Mass splittings in the Bg spectrum, converted to MeV using 1/at from Table 11. The bar 
over Bs indicates the spin-averaged mass {Mb, + )/4 was used. 


P 

l/ag (GeV) 

action 

K 

B*g - Bg 

B/n - Bg 

Bsi - B/q 

B/n - Bg 

Bgi — Bg 

1.8 

0.7 

1-link 

5 

34.9(10.2) 

442(56) 

12.6(4.7) 

416(54) 

430(57) 

1.8 

0.7 

AsqTad-tn 

4 

31.9(2.4) 

285(78) 

10.1(3.4) 

261(80) 

272(77) 

1.719 

0.8 

1-link 

3 

21.1(2.7) 

471(25) 

0.8(2.8) 

456(25) 

456(26) 

2.131 

1.0 

1-link 

4 

30.7(3.5) 

315(105) 

23.1(24.6) 

292(107) 

321(130) 

2.131 

1.0 

AsqTad 

4 

25.0(4.8) 

523(94) 

35.0(36.1) 

504(92) 

545(101) 

2.4 

1.2 

1-link 

6 

25.6(12.1) 

425(60) 

-9.4(22.2) 

406(57) 

398(63) 

2.4 

1.2 

AsqTad-tn 

6 

32.4(8.0) 

403(56) 

14.2(21.7) 

380(60) 

392(66) 
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Fit result 


1.8 



FIG. 2: Results of several fits to Propagator of Fig. 1 plotted vs. number of exponentials in the fit. Parameters 
shown are the non-oscillating ground state energy {aEo), the oscillating ground state energy {aEi), and the 
non-oscillating first excited state energy {aE 2 = aEo + aAE 2 ). 
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B propagator, AsqTad 
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FIG. 3: Example of the poor fits obtained 
mproved staggered quark, caused by the 
^extb The fit shown has y?,,„/DoF = 8.9. 




aE 



FIG. 4: Dispersion relation for free massless fermions. The clotted line shows the continuum dispersion 
relation , the dashed line shows the dispersion relation for the naive fermions sinh^ aif = sin^ ap, 

and the solid lines show the real part of the dispersion relation for the Naik action. Note that the solution 
of the Naik dispersion relation which most closely follows the continuum dispersion relation is purely real 
until the branch point near ap ~ 1. 
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FIG. 5: B meson propagator on the anisotropic 8^ x 48 lattice, where = 3.7 GeV, with a am = 0.04 
improved staggered quark and Mq = 7.0 nonrelativistic heavy quark. The lattice artifacts due to the 
temporal Naik term do not contaminate the fit. The 4-exponential fit plotted has x^ug/DoF = 0.87. 
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FIG. 6: Pion speed-of-light squared vs. momenta on the coarse (8^ x 20) lattice using several actions. The 
1-link (diamonds) and AsqTad (circles) results are ours, compared to the clover (crosses) and D234 (squares) 
results of [16]. 
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FIG. 7: Kinetic mass for the Bs meson on the 12^ x 32 lattice with the improved staggered action {am = 0.10) 
and the 1/M NRQCD action {aMo = 5.0). Computed using (58). The dashed line marks the experimental 
measurement Mb^ =5.37 GeV, and the solid lines show the range given perturbatively. 
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FIG. 8: Kinetic mass for the B* meson on the 12^ x 32 lattice with the improved staggered action {am = 0.10) 
and the 1/M NRQCD action {aMo = 5.0). Computed using (58). The dashed line marks the experimental 
measurement Mb* = 5.42 GeV, and the solid lines show the range given perturbatively. 
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FIG. 9: The hyperfine splitting between Bg and B* computed on quenched lattices. Our NRQCD-staggered 
results (from the isotropic lattices) are the circles. The squares come from [29], the fancy square from [27], 
and the diamond from [28] which all used an NRQCD-clover action, and the fancy cross comes from an 
NRQCD-D234 calculation [30]. All tune the heavy quark mass as in this work (see text for elaboration). 
For comparison, the experimental measurement is 47.0 ± 2.6 MeV [33]. Error bars are statistical only. 
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